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Abstract 

This describes a C version of a infeasible primal-dual algorithm for positive definite 
quadratic programming with box constraints, proposed by Voglis and Lagaris. We also 
give a straightforward .C() interface for R. 
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Note: This is a working paper which will be expanded/updated frequently. All suggestions 
for improvement are welcome. The directory gifi.stat.ucla.edu/boxqp has a pdf version, the 
complete Rmd file with all code chunks, the R and C code, and the bib file. 


1 Introduction 

The problem we address in this paper is minimization of 

f(x) = ^ x'Cx + d'x (1) 

over a < x < b (elementwise). Here C is positive definite, so / is strictly convex. Also we 
assume without loss of generality that a < b. As a special case we can have either oq = — oo 
or bj = -poo or both. 

The Lagrangian is 


1 





£(x, A, /i) 


( 2 ) 


= -x'Cx + d'x — X'(x — a) — //(& — x), 
and the Karush-Kuhn-Tucker (KKT) necessary and sufficient conditions for a minimum are 


Cx T d — ATyU — 0, 

( 3 ) 

A > 0, 

( 4 ) 

h > o, 

( 5 ) 

o' 

1 

( 6 ) 

n'(b — x) — 0, 

( 7 ) 

a < x < b. 

( 8 ) 


2 Algorithm 

A straightforward primal-dual active set algorithm to solve the KKT conditions was suggested, 
implemented, and tested by Voglis and Lagaris (2004). There is no proof that the algorithm 
ends in a final number of steps and is not subjected to cycling. Intermediate iterates will 
be infeasible, and there is no guarantee that the objective function decreases in each step. 
Nevertheless the method seems to work well in practice. In future work we will implement 
the more reliable, but much more complicated, feasible active set algorithm of Hungerlander 
and Rendl (2015). 

Suppose in iteration k we have (x^ k \ X^ k \ /i^). Define the index sets 

1. Step 1: 

L ^ = {i | xf^ < a,i or {xf^ = and A ^ >0)} 

U^ = {i | xf^ > bi or ( x[ k) = 6* and nf^ >0)} 

S ^ = {i | Oj < xf^ < bi or (xf^ = a* and A< 0) or (x^ = and < 0)} 

2. Step 2: 

x[ k+1 ^ = a,i and l^\ k+1) = 0 for all i G lS k \ 

xf +1) = bi and Af +1) = 0 for all i G U (k \ 

A? +1) = 0 and /jf +1) = 0 for all i G S (k) . 

3. Step 3: Solve 

Ax (k+1) + b = A (fc+1) - p, (fc+1) (9) 
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for the unknowns 


xf +1) with i G R (fc) , 
pf +1) with i G LJ (k \ 
Af +1) with i G L (fc) . 


4. Step 4: If (a;( fc+1 ), \( k+1 \ ^( fc+1 )) satisfies the KKT conditions, stop. Otherwise set k 
equal to k + 1 and go to step 1. 


3 Implementation 

Computation is done in C, using conventions which make it easy to call the programs from R, 
as well as calling them from C routines that have no access to R. So all arguments are passed 
by reference, and functions do not return any values. Both I/O and memory allocation are 
the reponsibility of the calling program, which can be either R or C. To move painlessly 
(well, not quite) from one-based to zero-based vector indexing in C we use the inline function 
VINDEXO, to move from column-major-one-based to row-major-zero-based matrices we use 
the inline function MINDEXO. I hope the compiler gets the inline message. 

Most of the computational effort goes into solving (9). As Voglis and Lagaris (2004) show, 
we find x 4 ' fc+1 ' ) with i G by solving a positive definite linear system of order equal to the 
number of elements in S^ k \ The computation of the dual variables A^’ +1 - ) and is then 

a simple matrix multiplication. 

To solve the positive definite linear system we use the LAPACK routine dposv(). There are 
many ways to get at dposv (). We can get it from the system LAPACK library, which means on 
OS X from the vecLib framework. Or we could use Intel’s mkl or similar. Instead we have cho¬ 
sen to install a separate LAPACK from HOMEBREW (2017) in /user/local/opt/lapack, 
and use the LAPACKE C-interface (LAPACKE (2016)). Your mileage may vary. 


4 Code 

4.1 R code 

dyn. load("boxqp. so") 

boxcqpR <- functionfa, b, lower, upper, itmax=10) { 
n <- nrow (a) 
x <- -solvefa, b) 
lambda <- rep(0, n) 
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mu <- rep(0, n) 
amat <- matrix (0, n, n) 
bvec <- rep (0, n) 
ind <- rep(0, n) 
h <- 
• C( 

"boxcqp" , 

as.double (a), 

as.double (b), 

x = as.double (x) , 

lambda = as.double (lambda), 

mu = as.double (mu), 

as.double (lower), 

as.double (upper), 

as.double (amat), 

as.double (bvec), 

as.integer (n), 

as.integer (ind), 

itel = as.integer (0), 

f = as.double (0) , 

as.integer (itmax) 

) 

return (list (x = h$x, f = h$f, itel = h$itel, lambda = h$lambda, mu = h$mu)) 

> 


4.2 C code 

#include <lapacke.h> 

#include <stdbool.h> 

inline int VINDEX(const int) ; 

inline int MINDEX(const int, const int, const int); 

inline int VINDEX(const int i) { return i - 1; } 

inline int MINDEX(const int i, const int j, const int n) { 

return (i - 1) + (j - 1) * n; 

} 

void boxcqp(const double *, const double *, double *, double *, double *, 

const double *, const double *, double *, double *, const int *, 

int *, int *, double *, const int *); 
void in(const double *, const double *, const double *, const double *, 
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const double *, const int *, int *, int *); 
void sb(double *x, double *lambda, double *mu, const double *iower, 
const double *upper, const int *ind, const int *n); 
void ss (const double *, const double *, double *, double *, double *, 

const double *, const double *, double *, double *, const int *, 
const int *, const int *); 

void co (const double *, const double *, const double *, const double *, 

const double *, const int *, const int *, bool *); 

void ff (const double *, const double *, const double *, const int *, double *); 

void dposv (const int *, const int *, double *, double *, int *); 

void mprint(const double *, const int *, const int *); 

void check (const double *, const double *, const double *, const double *, 
const double *, const double *, const double *, const int *, 
const int *); 

void dposv (const int *n, const int *nrhs, double *a, double *b, int ^status) { 
lapack_int info =0, nn = (lapack_int)*n, nr = (lapack_int)*nrhs; 
info = LAPACKE_dposv(LAPACK_COL_MAJOR, 1 U 1 , nn, nr, a, nn, b, nn); 

^status = info; 
return; 

> 

void boxcqp( const double *a, const double *b, double *x, double *lambda, 

double *mu, const double *lower, const double *upper, double *amat, 
double *bvec, const int *n, int *ind, int *itel, double *f, 
const int *itmax) { 
int nzero = 0, ktel = 1, isit = 0; 
bool *opt = (bool *)&isit; 
while (true) { 

(void)in(x, lambda, mu, lower, upper, n, ind, &nzero); 

(void)sb(x, lambda, mu, lower, upper, n, ind); 

(void)ss(a, b, x, lambda, mu, lower, upper, amat, bvec, n, ind, &nzero); 
(void)co(x, lambda, mu, lower, upper, n, ind, opt); 
if ((ktel == *itmax) || (*opt == true)) { 
break; 

} 

ktel++; 

> 

*itel = ktel; 

(void)ff(a, b, x, n, f); 
return; 

> 

void in(const double *x, const double ^lambda, const double *mu, 
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const double *lower, const double *upper, const int *n, int *ind, 
int *nzero) { 
int nvar = *n, nz = 0; 
for (int i = 1; i <= nvar; i++) { 
int ii = VINDEX(i); 
if (x[ii] < lower [ii]) { 
ind[ii] = -1; 

} 

if (x[ii] > upper [ii]) { 
ind[ii] = 1; 

} 

if ((x[ii] > lower [ii]) && (x[ii] < upper [ii])) { 
ind[ii] = 0; 
nz++; 

} 

if (x[ii] == lower [ii]) { 
if (lambda[ii] >= 0) { 
ind[ii] = -1; 

} else { 

ind[ii] = 0; 
nz++; 

} 

> 

if (x[ii] == upper [ii]) { 
if (mu [ii] >= 0) { 
ind[ii] = 1; 

} else { 

ind[ii] = 0; 
nz++; 

} 

} 

> 

*nzero = nz; 
return; 

> 

void sb (double *x, double *lambda, double *mu, const double *iower, 
const double *upper, const int *n, const int *ind) { 
int nvar = *n; 

for (int i = 1; i <= nvar; i++) { 
int ii = VINDEX(i); 
if (ind[ii] == -1) { 
x [ii] = lower [ii] ; 
mu[ii] = 0.0; 
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} 

if (ind[ii] == 1) { 
x[ii] = upper [ii] ; 
lambda[ii] = 0.0; 

} 

if (ind[ii] == 0) { 
lambda[ii] = 0.0; 
mu[ii] = 0.0; 

} 

> 

return; 

> 

void ss (const double *a, const double *b, double *x, double *lambda, double *mu, 
const double *lower, const double *upper, double *amat, double *bvec, 
const int *n, const int *ind, const int *nzero) { 
int nvar = *n, nz = *nzero, status =0, k = 1, nr= 1; 
for (int i = 1; i <= nvar; i++) { 
int ii = VINDEX(i); 
if (ind[ii] == 0) { 

int kk = VINDEX(k); 
bvec [kk] = -b[ii] ; 
int 1=1; 

for (int j =1; j <= nvar; j++) { 
int jj = VINDEX(j); 
double aa = a[MINDEX(i, j, nvar)]; 
if (ind[jj] == -1) { 

bvec[kk] -= aa * lower[jj]; 

> 

if (ind[jj] == 1) { 

bvec[kk] -= aa * upper [jj]; 

> 

if (ind[jj] == 0) { 

amat[MINDEX(k, 1, nz)] = aa; 

1 ++; 

> 

} 

k++; 

} 

> 

(void)dposv(nzero, &nr, amat, bvec, &status); 
k = 0; 

for (int i = 1; i <= nvar; i++) { 
int ii = VINDEX(i); 
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if (ind[ii] == 0) { 
x[ii] = bvec [k]; 
k++; 

} 

> 

for (int i = 1; i <= nvar; i++) { 
int ii = VINDEX(i); 
if (ind[ii] == -1) { 
lambda [ii] = b[ii]; 
for (int j =1; j <= nvar; j++) { 
int jj = VINDEX(j); 
double aa = a[MINDEX(i, j, nvar)]; 
lambda[ii] += aa * x[jj] ; 

} 

} 

if (ind[ii] == 1) { 
mu[ii] = -b [ii] ; 

for (int j =1; j <= nvar; j++) { 
int jj = VINDEX(j); 
double aa = a[MINDEX(i, j, nvar)]; 
mu[ii] -= aa * x [ j j ] ; 

} 

} 

> 

return; 

> 

void co (const double *x, const double *lambda, const double *mu, 

const double *lower, const double *upper, const int *n, const int *ind, 
bool *opt) { 
int nvar = *n; 
bool isit = true; 
for (int i = 1; i <= nvar; i++) { 
int ii = VINDEX(i); 
if (ind[ii] == 0) { 

if (x[ii] < lower [ii]) { 
isit = false; 

} 

if (x[ii] > upper [ii]) { 
isit = false; 

} 

} 

if (ind[ii] == 1) { 
if (mu [ii] < 0) { 
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isit = false; 

} 

} 

if (ind[ii] == -1) { 

if (lambda[ii] < 0) { 
isit = false; 

} 

} 

> 

*opt = isit; 
return; 

> 

void ff (const double *a, const double *b, const double *x, const int *n, 
double *f) { 
int nvar = *n; 
double sa = 0.0, sb = 0.0; 
for (int i = 1; i <= nvar; i++) { 
int ii = VINDEX(i); 
sb += b [ii] * x [ii] ; 
for (int j = 1; j <= nvar; j++) { 
int jj = VINDEX(j); 
double aa = a[MINDEX(i, j, nvar)]; 
sa += aa * x[ii] * x [ j j ] ; 

} 

> 

*f = sa / 2.0 + sb; 

> 
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